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RAYLEIGH-BENARD SIMULATION USING GAS-KINETIC BGK SCHEME IN THE 

INCOMPRESSIBLE LIMIT * 

KUN XUt AND SHIU-HONG LUI* 


Abstract. In this paper, a gas- kinetic BGK model is constructed for the Raylcigh-Benard thermal 
convection in the incompressible flow limit, where the flow field and temperature field are described by two 
coupled BGK models. Since the collision times and pseudo-temperature in the corresponding BGK models 
can be different, the Prandtl number can be changed to any value instead of a fixed Pr = 1 in the original 
BGK model. The 2D Raylcigh-Benard thermal convection is studied and numerical results are compared 
with theoretical ones as well as other simulation results. 

Key words, thermal instability, incompressible flow, gas-kinetic scheme 

Subject classification. Applied Numerical Mathematics 

1. Introduction. The use of a code for compressible flow to study incompressible fluid has attracted 
much attention in the past years. Since compressibility is proportional to the Mach number squared, Sp/p ~ 
M 2 , it is negligible once the Mach number is lower than 0.15. In many numerical test cases, such as the cavity 
flow, the results from compressible codes are almost identical to the results from incompressible codes[3, 9, 13] . 
It is also realized that using a compressible code for incompressible simulations have advantages. For example, 
a Poisson solver is avoided and parallelization of the code can be easily implemented. 

If thermal effects arc involved in the incompressible flow, a simple adaptation of a compressible code 
here bears potential danger. The reason is that the density varies with the temperature; this variation 
cannot in general be neglected, and therefore, even at small velocities, the density of a non-uniformly heated 
fluid cannot be supposed constant. For example, across the thermal boundary layer, the pressure is almost 
constant. If the temperature changes substantially, say by 10%, in the layer, then the energy equation will 
cause a 10% density change due to the ideal equation of state p = pRT. In reality, the density change 
is minimal with any reasonable temperature variation in the liquid. So, the compressibility effect is more 
severe in the thermal problem than that for the pure Mach compression problem where 5p/p ~ M 2 . It 
is certainly true that we can use other equations of state to describe a slightly compressible liquid. See 
[10] and references therein. There, the ability to recover the correct thermal effects is still questionable. 
In most current literature about the application of compressible codes to incompressible flows, thermal 
compressibility seems to be ignored. 

In order to reduce the compressibility in the compressible code for the thermal problem, we have to, 
in some ways, decouple the mass and momentum from the energy equation. In this paper, two pseudo- 
temperatures arc used to model the Raylcigh-Benard thermal convection problem in the incompressible 
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limit. In the current model, the velocity field and temperature field are described by two BGK models with 
different collision times. As a consequence, the Prandtl number can be changed to any value by modifying 
the collision times. 


2. Gas-Kinetic BGK Models for Rayleigh Benard Thermal Convection. In this section, we arc 
going to construct BGK models to study the following incompressible Navier-Stokes equations with thermal 
effect, 


( 2 . 1 ) 


r f£ + V-(pU) = 0, 

< f + U-VU = -5 +i/V 2 U-G, 
k + V • (pTU) = V • (kVT), 


where p is the density which is a constant in the incompressible limit, U the velocity, p the pressure, k the 
coefficient of thermal conductivity, and T the temperature. Note that pT is the thermal energy. For the 
Rayleigh-Benard convection in a two-dimensional box, the Boussinesq approximation gives 


pG - ppG 0 (T - T m )y, 

where G 0 is the gravitational constant, T m the average value of the top and bottom temperatures, y the unit 
vector in the vertical direction, and 0 the coefficient of volume expansion. For authoritative treatments of 
this problem, sec, for example, [2] and [6]. 

In order to the recover the above equations, gas-kinetic models can be constructed in the following forms, 

df f eq - f 

(2.2) ijr + u • V/ = — J -+F, 

Ol T„ 


(2.3) 


dh h eq -h 

— + u ■ V/i = , 

Ot Tr 


where u = (u, v) is the x and y components of the particle velocity. Eq.(2.2) is used to recover the mass 
and momentum equations, and also the velocity flow field. Eq.(2.3) is for the thermal energy evolution. The 
equilibrium states f eq and h eq have the following forms 


f e <! =zp(- i) e -M( u “ u ) 2 ), 
7 r 


h eq = pT{ — )e- A2 « u - u > 2 > 

7 r 

where Ai and A 2 can be expressed as 


Ai = 


1 

2 RTi 


and 


A2 = 


1 

2 RTv 


with the two pseudo- temperatures T\ and T 2 . Here T is the real temperature to be simulated. Note that 
T\ and T 2 arc both constants in the current model, and the value of either T\ or Ai determines the artificial 
sound speed of the flow field. In the above BGK models, the compressibility is determined from Eq.(2.2) 
with the equation of state p = pRT\, which is totally decoupled from the real temperature T. The external 
forcing term F in Eq.(2.2) can be approximated as [7], 


F — 2AG ■ (u — U)/ e<7 , 


TT 
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from which the buoyancy force can be recovered. 

In the course of particle collisions, the compatibility condition is satisfied in the BGK models, 


I ( r - /) 


/i\ 

u 

w 


dudv — 0, 


and 


J ( h eq — h)dudv — 0. 


By using Chapman-Enskog expansion, Eq.(2.1) can be recovered exactly in the incompressible limit, 
with the viscosity coefficient 


v = r v pRT\ 


and the heat conduction coefficient 


k = T c pRT 2 . 


Different from the original BGK model [1], here both coefficients are decoupled from the fluid temperature 
T. As a result, the Prandtl number Pr becomes 



Ty_ T\_ 

t c T 2 ’ 


which can be changed to any value by choosing different r v , r c , T \ , or T 2 . 


3. Numerical Scheme for the BGK Models. For a finite volume scheme, we need to evaluate the 
numerical fluxes across a cell interface, and the flux function depends on the gas distribution function. In 
this section, the BGK scheme to solve Eq.(2.2) and (2.3) for fluxes will be presented. 

Firstly, for Eq.(2.2) we arc going to use the operator splitting method to solve the equation into two 
steps 

feq _ f 

(3.1) ft + Uf x + Vfy — , 

T v 


and 


(3.2) 


ft=F. 


For Eq.(3.1), in the smooth incompressible limit, the general solution of / in the above equation at the cell 
interface and time t can be simplified as [15], 


(3-3) 


1 f* 

f{*i+i/ 2 ,j,Vi+i/ 2 ,j,t,u,v) = — / f eq (x',y',t',u,v)e~^ t )lT "dt' 

T v J-o o 


where x* = £*+1/2 t j — u(t — t r ) and y f = y x +\f2,j — v{t — t f ) is the trajectory of a particle motion. Generally, 
the equilibrium state f eg around the center of the cell interface (#14-1/2, j = #o ? Vi+\/2,j — yo) 5 and the initial 
time step (t = 0) can be approximated as 


(3.4) 


f eq (x, y, t,u, v) = (1 + (x - z 0 )a + (y - y 0 )b + tA) g Q) 
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where go is the local Maxwellian located at the center of a cell interface, 

(3.5) g 0 = Po (^-) e -M(«-t/o) 2 +(«-v„) 2 ] 

Note again Ai is a constant. The dependence of a, b,A in Eq.(3.4) on the particle velocities can be obtained 
from the Taylor expansion of a Maxwellian and have the forms 


a — a i 4- a 2 u + a 3 u 

= (-^ + 2A 1 f/ 0 ^ + 2AiV 0 ?) - ZMUojf-u - 2X 1 V 0 ^-v, 

po OX OX OX OX OX 

b = 6i + b 2 u -I- b$v 

= (— + 2X 1 Uo^- + 2A!V 0 ^) - 2X 1 U 0 ^u - 2A 1 V 0 ^i>, 

Po dy oy dy dy dy 

A = A\ 4- A 2 u -f A 3 v 

= { Jo% + 2XiUo w + 2XiVo % ] - 2XlUo l£ U ~ 2XlVo W V ’ 

where all parameters (dp/dx, dU/dx, dV/dx) and (dp/dy,dU jdy,dV/dy) at t = 0 can be obtained from 
the initial reconstructions of the macroscopic variables dp/dx, dpjdy, d(pU)/dx,.,, For example, a 2nd-order 
interpolation gives 


Po — ~ (Pij + Pi+ l,j) 


Uo = + (pU)i+i t j) 

Vo = — ((pV)iJ + (pV)i+lj) 


dp 

dx 


1 

Ax 


(p*+ l,j pi,j ) 


dp 

dy 


+ Pi,j+ 1) “ 2 (P'+hj-i + 


where Ax, Ay are the cell sizes in the x and y directions. 

After substituting Eq.(3.4) into Eq.(3.3), the final gas distribution function at a cell interface is 

(3.6) f(x 0 ,yo,t,u,v) = g 0 ( 1 - r v (ua + vb) + (t -t„)A). 

The only unknown in the above equation is A , which depends on dp/dt,dU/dt and dV/dt. Since 

f eq (x Ql y 0 , t , u, v) = go(l + At), 
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together with the compatibility condition 


u dudv — 0, 


along time t and at x = A can be uniquely determined from 


J g 0 (ua + vb + A) 


u dudv = 0, 


r dt 

JL m 

n dt 

Po l 

\ dt 


\ w f 1 ) 

) po w 


( a\ < u > +a 2 < u 2 > +03 < uv > +b\ < v > +62 < uv > +63 < u 2 > \ 

= - a\ < u 2 > +a 2 < u 3 > -f a 3 < u 2 u > +&i < uu > +& 2 < uv 2 > +&3 < uu 2 > > 

< uv > -fa 2 < > +^3 < uv 2 > +&i < v 2 > +62 < uu 2 > +63 < v 3 > J 

where the detail formulation of < u n u m > can be found in the Appendix. Therefore, the above equation 
uniquely determines dp/dt , dU/dt and dV/dt , so A is obtained. 

After determining / in Eq.(3.6), the time-dependent numerical fluxes in the x-direction across the cell 
interface can be computed as 

/m t n 

(3.7) T p u = I u\ u g 0 {l + r l/ {au+ bv) + (t - T v )A)dudv. 

K T PV ) i+l/2,} \ V ) 

Once again, the moments of u and v can be easily obtained from the recursive relation shown in the Appendix. 
By integrating the above equation for a time step At, we get the total mass, momentum transport. Similarly, 
Gij+ 1 / 2 ? the fluxes in the y direction can be obtained by repeating the above process in the y direction. With 
both fluxes in the x and y directions, we can update the flow variables inside each cell (i, j) by 


= PU 

\pV 


r* 1 1 1 

J ° V \p n pGo{T n - Tm) 


where the effect from Eq.(3.2) has been accounted for in the above equation. 

Once Eq.(2.2) is solved, the scheme for Eq.(2.3) can be constructed similarly. For example, we can 
expand h eq as 


h eq {x, y , t, u, v) = fr 0 (l + (x - x o) a h 4 (y - yo)bh, + tA h ), 


where 


at a cell interface, and 


fro = (P 0 T 0 ) 


-\i{{u-U 0 ) 2 + {v-V 0 f) 


a h = a h\ + Q'h2'U + &h3V 

= ( 1 d [ pT) + 2A 2 U 0 ^~ + 2A 2 V 0 jb - 2A 2 U 0 ^u - 2X 2 V 0 ^v, 

p 0 T 0 dx dx ox ox ox 
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bh = bhi + bhiu + bh 

1 d( P T) 


= ( 


P 0 T 0 


- f 2A2C/0 t; — h 2A2Vo-^— ) — 2 X 2 U 0 —— u — 2A2V0 -jj— u, 

5?/ dy dy dy dy 


Ah, — Ah,\ + 


f 1 flOT 

poTo <9t 


+ 2A 2 C/ 0 ^ + 2 ^V 0 ^) - 2 *2Uo~u - 2X 2 V 0 ^v, 


which arc closely related to the coefhcients of a, b and A. In other words, the evolution of h is not totally 
independent of the evolution /, and dU/dx, dV/dx, ... in the above equations are the same as the corre- 
sponding terms in the equations defining a, b t A earlier. Hence, the only unknowns are T 0 ,dT/dx,dT /dy 
and dT/dt. In order to determine all unknowns, at t = 0, the following interpolations can be used to get 
p 0 T 0 and d(pT) / dx , d(pT) /dy. The linear reconstruction of thermal energy pT is necessary with 


PoTo — 0.5((pT)i ? j -b (pT)i+ij), 


and 


d(pT) 

dx 


— (( P T), +lti - (pT)ij), 


d(pT) 

dy 


2^((A)^b)»+i/2,i+i - (PoT 0 )i + i/ 2 j-i). 


The final solution of h at the center of the cell interface is 


(3.8) 


h(xo,yo,t,u,v) = h 0 (l - + u&h) + (t - r c )A h ), 


and the dT/dt term in Ah is determined by applying the compatibility condition 

J ( h eq — h)dudv = 0, 


along (xo,yo,t) y which similarly gives 

J Ahhodudv = — J {ahu + bhv)hodudv. 

Once h is determined in Eq.(3.8), the numerical flux for the thermal energy is 


T pT = 


J uhdudv , 


and the thermal energy inside each cell can be subsequently updated. 

4. Results. The Rayleigh- Benard problem offers a first approach to a complicated convective flow. In 
this case, with the gravitational force in the vertical direction a horizontal layer of viscous fluid is heated from 
the bottom while the top boundary is maintained at a lower temperature. When the temperature difference 
between the top and bottom boundaries is increased above a certain threshold, the static conduction state 
becomes unstable to any small disturbance and the system become convective. 

In our calculations, the horizontal and vertical length scales are L = 2.0 and H = 1.0, respectively. The 
temperatures at the bottom and top are Tbottom = 1-0 ,T top = 0.0, with the difference AT = 1.0. Non-slip 
boundary conditions are implemented at the bottom and top boundaries by reversing the flow velocities in 
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the ’ghost’ cell next to the simulation domain. Periodic boundary conditions are used for the temperature 
along the sides of the box. In our current study, we fix Go = 1.0 and (3 = 0.1. 

The Rayleigh number is defined as 

n 0ATG o H 3 
R ~ ^ • 

From the above relation and Pr = i s/k, the viscosity coefficient can be determined: 


/3ATG 0 H 3 Pr 

u ~ V R 

Consequently, the collision time r v in Eq.(2.2) is fixed with 


r v — 2 Ail/, 

and r c in Eq.(2.3) is 

Tc “ 

AiPr 

Since in the simulations, the CFL time step At is almost a constant, in order to keep the collision time r v to 
be around 10“ 1 At, we have to choose Ai properly. In most calculations, Ai is on the order of 10 _1 . Although 
the numerical scheme is general for any Pr, we used P r — 1, Ai = A 2 and r v = r c . 

As a first test, we tried to get the critical Rayleigh number for the onset of thermal convection. With a 
80 x 40 mesh, we have simulated this problem with two supercritical Rayleigh numbers R — 1720 and R — 
1735 separately. In each case, we calculate the maximum y-component velocity in the whole computational 
domain at each time step. The time- dependent amplitude of the y- velocity on a 80 x 40 mesh is shown in 
Figure 5.1, from which we can estimate the critical Rayleigh number by fitting the curve to V ~ cxp(a(P — 
R c )t) y where R c is the critical Rayleigh number. From the exponential growth rates, we found that the 
critical Rayleigh number in our calculations is R c = 1711.17, which is 0.22% away from the theoretical value 
1707.76 (which is actually for a box of width 2.0158). For other meshes, the calculated critical Rayleigh 
numbers are listed in Table 1. 

Table 1. Critical Rayleigh numbers calculated on different meshes. The error is calculated relative to the 

theoretical value. 


Grid Size 

Ra c 

Error 

20 x 10 

1756.22 

2.84% 

40 x 20 

1729.43 

1.27% 

80 x 40 

1711.45 

0.22% 

theory 

1707.76 



Once the Rayleigh- Benard convection is stabilized, the heat transfer between the top and bottom is 
greatly enhanced. The enhancement of the heat transfer can be described by the Nussclt number, 

, < VT > 

u_1+ kAT/H' 

where V is the vertical velocity, AT is the temperature difference between the bottom and top walls, H 
is the height of the box, and < ... > represents the average over the whole flow domain. Figure 5.2 is the 
calculated relationship between the Nusselt number and the Rayleigh number. The simulation results by 
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Clever and Bussc [4] is also included. As shown in the figure, our results are very close to those by Clever 
and Busse. But, at higher Rayleigh numbers, our values of the Nusselt number is a little bit smaller than 
that in [4], and thus underestimating the amount of heat transfer. Similar results are obtained using lattice 
Boltzmann methods [12, 8], 

Typical temperature and stream function contours arc shown in Figures 3-8 with Ra = 5, 000, 10,000 
and 50,000. As the Rayleigh number increases, two trends were observed for the temperature distribution: 
enhanced mixing of the hot and cold fluids and an increase in the temperature gradients near the bottom 
and top boundaries. Both trends enhance the heat transfer in the box. 

As another benchmark problem, we have tried one case in [5]. This problem is that of the two-dimensional 
Boussinesq flow in a square with H = L = 1.0 and Prandtl number Pr ~ 0.71, which is done by setting 
Ai = A 2 and r c = r u jPr in our code. Both velocity components arc zero on the boundaries. The horizontal 
walls are insulated, and the vertical sides are at temperatures = 1.0 and T r i g ht — 0.0. In this case, the 
Nusselt number is defined as 

1 < UT > 

N *- 1+ kAT/L ' 

The results for the streamline and temperature contours at R — 10 5 are shown in Figure 5.9 and 5.10. 
With R = 10 5 , the average Nusselt number in the whole domain is listed in Table 2 for different mesh 
sizes. Contrary to the last test case, our result overestimates the heat transfer. A larger Nusselt number is 
obtained. 

Table 2. Nusselt numbers calculated on different meshes. The error is calculated relative to the numerical 

result in [5]. 


Grid Size 

Nusselt Number 

Error 

20 x 20 

4.590 

1 . 77 % 

40 x 40 

4.563 

1.17% 

80 x 80 

4.540 

0.66% 

reference [5] 

4.510 



5. Conclusion. In this paper, a two-temperature gas-kinetic BGK model for convective thermal flow 
is constructed. A numerical scheme has subsequently been developed. As an application, the 2D Rayleigh- 
Benard case is studied. The simulation results arc very close to those obtained by other methods. To study 
the incompressible flow phenomena using the compressible model is an attractive research area. In order to 
simulate thermal effects in an incompressible fluid, the decoupling of the energy equation from the mass and 
momentum equations seems necessary, because the relation between temperature and volume changes are 
different for incompressible and compressible fluids. Compared with the lattice BGK methods, the current 
approach with continuous particle velocity has advantages in terms of stability and efficiency. The time step 
used in the current method is the CFL time step which is about one order of magnitude larger than the 
particle collision time which is usually used in the lattice BGK method[ll]. 

In this paper, the temperature evolution equation only includes advection and diffusion terms. The 
viscous heating term in the Navier-Stokes energy equation is ignored due to the simplicity of the model. 
The construction of a two-temperature BGK model with the viscous heating term in the thermal energy 
evolution equation is an interesting and important problem. The research in this direction will help us to 
find an efficient kinetic scheme to simulate incompressible flow, and pave the way to simulate a flow mixing 
compressible gas and incompressible liquid. 




Acknowledgments. We would like to thank Dr. X. He for his helpful discussion and advice in this 
work, and thank Dr. L. Luo for his valuable comments. 

Appendix: Moments of the Maxwellian Distribution Function. In the gas-kinetic scheme, we 
need to evaluate moments of the Maxwellian distribution function with unbounded integration limits. Here, 
we list some general formulas. 

Firstly, we assume that the Maxwellian distribution for a two-dimensional flow is 

g = p ^) e -H(u-U)*+(v-Vf)' 

7 r 

Then, by introducing the following notation for the moments of g , 

p < ... >= J (...)gdudv, 

the general moment formula becomes 


< u n v m >=< u n >< V th >, 

where n, m are integers. When the integration limits arc from — oc to +oo, we have 

< u° >= I 

< u >= U 


< u"+ 2 >= TJ < u n+1 > 1 < u n > . 

lA 

Similarly, 


< v° >= 1 

< v >= V 


< V m+2 >= V < v m+1 > +^r-i <v m > . 
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Fig. 5.1. Time history of the maximum vertical velocities. 


& 4h 
E 


+: Currant BGK Scheme 
O: Clever & Busse 


I 1 - 

io 3 


10 

Rayleigh Number 
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